The structure and phase stability of CO adsorbates on Rh(llO) 
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The structure of CO adsorbates on the Rh(110) surface is studied at full coverage using first- 
principles techniques. The relative energies of different adsorbate geometries are determined 
by means of accurate structure optimizations. In agreement with experiments, we find that a 
p2mg(2 x l)-2CO structure is the most stable. The CO molecules sit on the short-bridge site (car- 
bon below) with the molecular axis slightly tilted off the surface normal, along the (001) direction. 
Configurations corresponding to different distributions of tilt angles are mapped onto an anisotropic 
2D Ising model whose parameters are extracted from our ab-initio calculations. We find that an 
order-disorder phase-transition occurs at a temperature T c ~ 350 °K. 
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I. INTRODUCTION II. COMPUTATIONAL METHOD 



Rhodium surfaces are attracting a wide scientific and 
technological interest due to their catalytic properties, 
particularly because they act so as to reduce the energy 
activation barrier for the reaction 2CO + 2NO — > 2CO2 + 
N2, and thus to eliminate the two poisonous CO and NO 
gases from the pollution emission of combustion engines. 

The stable structure of the Rh(110) clean surface is 
unreconstructed. However, if prepared in a convenient 
way with oxygen adsorption it may also present — upon 
thermal desorption — meta-stable (1 X n), (n — 2, 3, 4, 5), 
missing- or added-row structures which revert to the un- 
reconstructed one at temperatures above 480 °K 

The adsorption of CO molecules on Rh(110) has been 
studied by means of a variety of techniques p|~|L2]] . The 
adsorption of 1 monolayer (ML) of carbon monoxide on 
the unreconstructed surface results in a (2 x l)p2mg 
structure, with the C atom bound in the short bridge 
sites along the (110) direction, and the molecular axis 
alternatively tilted with respect to the surface normal, 
towards (001) direction ||. In Ref. [§, the (2 x l)p2mg 
LEED pattern was reported to disappear at temperatures 
higher than » 270-^280 °K, well below the desorption of 
CO from the surface. This fact was tentatively explained 
in terms of an order-disorder phase-transition. 

In this paper the structure and phase stability of one 
ML of CO molecules adsorbed on the Rh(110) (1 x 1) 
surface are studied from first principles and by mapping 
the low-lying energy configurations corresponding to the 
different distributions of tilt angles onto an anisotropic 
2D Ising model. The latter is then simulated using a 
standard Metropolis Monte Carlo algorithm. The order- 
disorder transition temperature estimated from our sim- 
ulations is 350 ± 50 °K , in fair agreement with exper- 
imental findings [||. 



Our calculations are based on density functional theory 
within the local-density approximation (LDA) |L3| . |l4| . 
using Ceperley- Alder correlation energies |]l5| . The one- 
particle Kohn-Sham equations are solved self-consistently 
using plane-wave (PW) basis sets in a pseudopotential 
scheme. Because of the well known hardness of the norm- 
conserving (NC) pseudopotentials for the O and — to a 
lesser extent — Rh atoms, we make use of ultra-soft (US) 
pseudopotentials ]l6| which allow an accurate descrip- 
tion of the O and Rh valence pseudo-wavefunctions with 
a modest basis set including PW's up to a kinetic-energy 
cutoff of 30 Ry. In the case of Rh, we found it conve- 
nient to treat the s and p channels using a NC poten- 
tial, while the US scheme is applied only to the hard 
d orbital |i~7 18 1. With such a small basis set, the ac- 
curacy is slightly improved if also C is treated within 
the US scheme, which we decided to do. Brillouin- 
zone (BZ) integrations are performed using the Gaussian- 
smearing |l9|] special-point [ p0| technique. In agreement 
with Ref. JL8|, we find that the structural properties of 
bulk rhodium are well converged using a first-order Gaus- 
sian smearing function Jl| of width g = 0.03 Ry and 10 
special k-points in the irreducible wedge of the BZ (IBZ). 
The isolated surface is modeled by a periodically repeated 
supercell. We have used the same supercell for both the 
clean and the CO-covered surfaces. For the clean surface 
we have used 7 atomic layers plus a vacuum region corre- 
sponding to w 9 layers. For the CO-covered surface the 
7 Rh layers are completed by one layer of CO molecules 
on each side of the slab: in this case the vacuum region is 
correspondingly reduced to ~ 5.5 atomic layers. We have 
used the same Gaussian-smearing function as in the bulk 
calculations with 8 special k-points in the surface IBZ. 
Convergence tests performed with a value of a twice as 
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small and a correspondingly finer mesh of special points 
resulted in no significant changes in total energies and 
equilibrium geometries. The latter are found by allowing 
all the atoms in the slab to relax until the forces acting 
on each of them are smaller than 0.5 x 10~ 3 Ry/ao- 

III. RESULTS 
A. Structural analysis 

The clean Rh(110) surface is unreconstructed. An 
analysis of LEED data suggests that the top interlayer 
spacing is reduced by 6.9± 1.0% relative to the bulk inter- 
layer spacing, while the second interlayer spacing would 
expand by 1.9 ± 1.0% Q. Our ab initio data indicate 
a relaxation of —9% and of 3.5% in the first and second 
interlayer spacings respectively. 

For the CO-covered surface, LEED data indicate that 
the molecules are bound in the short bridge site between 
two first layer rhodium atoms in the (001) direction with 
the molecular axis tilted by 24 degrees from the surface 
normal, forming a (2 x l)p2mg structure. In princi- 
ples there are many different ways to arrange the CO 
molecules so as to obtain the same LEED pattern. We 
concentrate our attention on three possible adsorption 
sites: i) the short-bridge one described above; ii) the on- 
top one, in which the CO molecule is located on top the 
first-layer atoms; and Hi) the hollow site, formed by two 
first-layer atoms in the (001) direction and one second- 
layer atom. In agreement with the outcome of the LEED 
analysis, we find that the short-bridge site is the most fa- 
vorable. The relative energies of the other two sites with 
respect to the short bridge — assuming a 1 x 1) structure 
in all cases — are: 0.19 eV (hollow) and 0.34 eV (on top). 
We find that the angle between the surface normal and 
the CO molecular axis is a — 17 ± 2 degrees, and that 
the angle between the Rh-C bond and the surface normal 
is S = 13 ± 2 degrees; the Rh-C bond length is 2.02 A 
and the C-0 distance is 1.17 A. The rhodium substrate 
presents an outward relaxation of the first layer of 2.8% 
with respect to the bulk interlayer spacing. These re- 
sults are summarized in Table |lj together with similar 
ones obtained for six other different surface geometries 
(See Fig. @). 

From Table § we see that the (lxl) and (1 x 2) ge- 
ometries are degenerate within our error bar which we 
estimate to be ±1 meV, and that the uncertainty about 
the corresponding tilt angle is very large. This behav- 
ior can be understood by a simple qualitative model of 
the surface energetics which also accounts for the ob- 
served ordering of the structures. In order to disentan- 
gle the relative importance of the adsorbate-substrate 
and adsorbate-adsorbate interactions, we have modeled 
the former by a (2 x 2) supercell in which a single CO 



TABLE I. Structural data for seven different surface struc- 
tures (see Fig. 1). a and S are respectively the angles between 
the surface normal and the C-0 and the Rh-C axis, AE is the 
energy difference per molecule between the (n x m) and the 
(2 x 1) structures. The theoretical error is estimated to be 
~ 2 meV. The experimental values refer to the (2 x 1) struc- 
ture 
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24±4 a 


13±4 a 


1.13±0.09 a 


1.97±0.09 a 




2x1 


17 ±2 


13 ± 2 


1.17 


2.02 


0.0 
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<io 


< 10 


1.17 


2.02 


33.5 


1 x 2 
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<5 


1.17 


2.02 


33.5 


2x2 


13 ±2 


11 ± 2 


1.17 


2.02 


13.5 


2x2' 


13 ±2 


11 ±2 


1.17 


2.02 


21.5 


4 x 1 


16 ±2 


13 ±2 


1.17 


2.02 


17.0 


4x1' 


16 ±2 


12 ±2 


1.17 


2.02 


17.0 



From Ref. § 



molecule is constrained to sit at the same short-bridge 
site which would be preferred at full coverage. 

We observe that the dependence of the adsorption en- 
ergy on a is very weak up to a ~ 10°, and that it be- 
comes very steep above this angle. When the coverage 
increases, the dipole-dipole interaction becomes impor- 
tant and accounts qualitatively for the energy ordering 
of the structures displayed in Fig. 1. In the (2 x 1) struc- 
ture, nearest-neighbor molecules are tilted by opposite 
angles around the axis joining them (the (110) direc- 
tion) , while next-nearest-neighbor molecules are tilted by 
a same angle about the (001) axis which joins them. The 
dipole-dipole interaction favors both these arrangements 
of angles. The (2 x 2) geometry is similar to the pre- 
vious one as regards the nearest-neighbor interactions, 
whereas it is unfavored regarding next-nearest-neighbor 
interactions. The next higher energies are those of the 
(4x1) and (4 x 1') structures which are almost degen- 
erate because they have the same number of unlike tilt 
angles along the (110) row. The next structure is the 
(2 x 2') one which is characterized by an alternating ar- 
rangements of energetically favored and disfavored rows 
and columns of CO molecules. Finally, in the (1 x 1) 
and (1x2) structures the nearest-neighbor molecules are 
tilted by a same angle and the corresponding dipolar in- 
teraction is therefore independent of cv; it is only the 
weaker next-nearest-neighbor interaction which depends 
on a, more so for the (1 x 2) structure for which the sign 
of the dipole-dipole interaction energy is the same as that 
of the adsorbate-substrate interaction, while the two in- 
teractions tend to cancel for the other structure. In both 
cases, this behavior results in a very weak dependence 
of the energy upon a, and in an energy degeneracy of 
the two structures, within our error bars. We have also 
calculated the adsorption energy of the CO molecules de- 
fined as E^ G ° - ~ where Ef^ co is the 

total energy of the CO covered surface, E^ b is the to- 
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2x2' 4x1 4x1' 

FIG. 1. The seven different surface structures referred in 
the text and the clean Rh(llO) surface. 

tal energy of the clean Rh(llO) surface, and Ef^ co 
is the total energy of the CO, all the calculations be- 
ing done using the same slab geometry and the same set 
of k-points. The calculated adsorption energy is of 2.78 
eV/molecule, which has to be compared to the experi- 
mental value 1.1 eV/molecule Q|. This large discrepancy 
is a common feature of the LDA which is well known to 
overestimate absolute binding energies, whereas equilib- 
rium geometries and energy differences among them are 
usually predicted with a much higher accuracy (of the 
order of a few percents). 



spirit of the cluster expansion of the energy landscape 
of an alloy Q , the energy of each tilt-angle configura- 
tion can be expressed in terms of polynomials in the cr's. 
Because, of symmetry, odd-power polynomials are ab- 
sent from the cluster expansion. Restricting ourselves to 
second-order polynomials (spin-pair interactions) and ne- 
glecting all the couplings beyond next-nearest-neighbors, 
the cluster expansion of the surface energy would read: 

i,j \ S=±l 

+ J V X! <Jl '3+ S + ^ a i+S,3+S' J • (1) 

<5=±1 S,8'=±l J 

It is straightforward to see that: 

£■2x1 = Jy + Jx — J2', Ei x i — J x + J y + J2; 
Elx2 — Jx ~ Jy — J2', E2x2 = J2 — Jy — Jx', 
E 4X 1 = i?4xl' = Jy] E 2 x2' = 0, (2) 

where the subscripts refer to the structures of Fig. |]. 

i?2xi is the ground-state energy which we take as the 
reference energy. The (4 x 1) and (4 x 1') structures are 
degenerate within the present model, and their energy 
difference provides therefore an estimate of longer-range 
or many-spin interactions which have been neglected. 
Out Eqs. (2), one can extract four independent energy 
differences, which are linear functions of the three pa- 
rameters J Xl J y , and J 2- By disregarding one of these 
equations in turn, one obtains 4 different linear sistems 
for the J's which provide estimates for these parameters, 
which coincide within 1.5 meV. The average of the 
four set of parameters so obtained is: J x — 13.6 meV, 
J y = —3.4 meV, and J2 = 3.4 meV. Note the large dif- 
ference between the absolute values of J x and J v , which 
is due to a stronger coupling in the 'ziz-zag' (110) direc- 
tion, where the distance between neighboring molecules 
is smaller by a factor \[2 than in the orthogonal direction. 



B. Finite-temperature properties 

1. Mapping onto a 2D Ising model 

From Table | we see that the energy necessary to tilt 
the angle of a molecule is of the order of 10 -j- 30 meV, 
whereas the energy difference between different adsorp- 
tion sites is typically ten times as large. This fact 
indicates that — given the adsorption sites of the CO 
molecules — a well defined energy would correspond to 
each distribution of tilt angles, where every molecule is 
labeled by a variable which indicates in which direction 
it is tilted (the sign of the tilt angle, which we associate 
to an Ising-like variable, a — ±1). Much in the same 



2. Monte Carlo simulations 

The thermal properties of our system are obtained 
by standard Metropolis Monte Carlo simulations of the 
above Ising model (see for example Ref. [^3)). To this 
end, we have used a 40 x 40 square lattice with peri- 
odic boundary conditions. The order parameter of the 
transition between the (2x1) ordered phase and the 
disorderd phase where the tilt angles are distributed at 
random, is the Fourier coefficient of the spin-spin correla- 
tion function, M(q) = -4 J2 r e * q ''(fr^o)) & t wavevector 
q = (7r,0). The order-disorder transition temperature, 
T Cl is estimated looking at the maximum of the specific 
heat C. We have not attempted any finite-size scaling, 
but we have verified that the location of the transition 
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FIG. 2. Specific heat and Fourier transform of the corre- 
lation length for the 2D spin model (see text). 
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temperature is rather insensitive to the choice of the size 
of the system, by making a few simulations for a 64 x 64 
spin matrix. In Fig. || we show the behavior of the specific 
heat, C, and the order parameter, M(jr, 0), as functions 
of temperature. We estimate the critical temperature to 
be T c ps 350 ± 50 °K. The estimate of the error bar is 
based on the uncertainties on the coupling coefficients of 
the Ising model. 
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